% Part 3: Variational Methods (3D-Var, 4D-Var, DRP-4DVar)
% This file is included by bayesian_vs_3dvar_cn_main.tex

\section{3D-Var 理论概述}
GSI 的 3D-Var 代价函数可表示为
\begin{equation}
J = \tfrac{1}{2}(x_a - x_b)^T B^{-1}(x_a - x_b) + \tfrac{1}{2}(H x_a - o_o)^T O^{-1}(H x_a - o_o) + J_c,
\end{equation}
其中 $x_a$ 为分析场，$x_b$ 为背景场，$B$ 为背景误差协方差，$H$ 为观测算子，$o_o$ 为观测向量，$O$ 为观测误差协方差，$J_c$ 为可选约束项。定义分析增量 $x = x_a - x_b$ 与创新向量 $o = o_o - H x_b$ 后，上式化为
\begin{equation}
J = \tfrac{1}{2} x^T B^{-1} x + \tfrac{1}{2}(H x - o)^T O^{-1}(H x - o) + J_c.
\end{equation}
为改善收敛性，引入控制变量变换 $y = B^{-1/2}x$，使得新变量在统计上无相关、方差为 1，有利于预条件共轭梯度 (PCG) 求解。GSI 通过平衡算子构建控制变量（如流函数、速度势、非平衡温度与地面气压等），并可混合静态与集合协方差获得流依赖误差信息。约束项 $J_c$ 用于实施数字滤波、湿度限制或弱约束等物理条件。

\section{四维变分资料同化 (4D-Var)}
三维变分(3D-Var)仅同化单一时刻观测，而四维变分(4D-Var)在时间窗口内同化多个时刻的观测，通过数值预报模式的切线性和伴随模式传递观测信息。

增量 4D-Var 的代价函数为：
\begin{equation}
J(\delta\mathbf{x}) = \frac{1}{2}\left[\delta\mathbf{x} - (\mathbf{x}^b - \mathbf{x}^g)\right]^T \mathbf{B}^{-1}\left[\delta\mathbf{x} - (\mathbf{x}^b - \mathbf{x}^g)\right] + \frac{1}{2}\sum_{i=0}^{n}\left[\mathbf{H}_i\delta\mathbf{x}(t_i) - \mathbf{d}_i\right]^T \mathbf{R}_i^{-1}\left[\mathbf{H}_i\delta\mathbf{x}(t_i) - \mathbf{d}_i\right],
\end{equation}
其中：
\begin{itemize}
\item $\mathbf{x}, \mathbf{x}^b, \mathbf{x}^g, \delta\mathbf{x} = \mathbf{x} - \mathbf{x}^g$ 分别为状态向量、背景场、初猜场、分析增量
\item $\mathbf{d}_i = \mathbf{y}_i^o - \mathbf{H}_i[\mathbf{x}^g(t_i)]$ 为 $t_i$ 时刻观测创新
\item $\mathbf{y}_i^o$ 为 $t_i$ 时刻观测，$\mathbf{H}_i$ 为观测算子，$\mathrm{H}_i$ 为切线性观测算子
\item $\mathbf{M}_i$ 为非线性预报模式，$\mathrm{M}_i$ 为切线性预报模式
\item $\mathbf{B}$ 为背景误差协方差矩阵，$\mathbf{R}_i$ 为 $t_i$ 时刻观测误差协方差矩阵
\end{itemize}

{\color{red}\textbf{[4D-Var 的时间维度]：与 3D-Var 仅在初始时刻约束状态不同，4D-Var 通过模式积分将初始增量 $\delta\mathbf{x}$ 传播至各观测时刻 $t_i$，利用整个时间窗口的观测约束初始条件。这需要构建切线性模式 $\mathrm{M}_i$ 和伴随模式 $\mathrm{M}_i^T$，计算梯度 $\nabla J = \sum_i \mathrm{M}_i^T \mathrm{H}_i^T \mathbf{R}_i^{-1}(\cdots)$。伴随模式的构建和维护是 4D-Var 最大的技术挑战。}}

4D-Var 相比 3D-Var 的优势：
\begin{itemize}
\item \textbf{时间一致性}：初始条件通过模式动力学约束，确保分析场满足物理平衡
\item \textbf{观测利用率}：同化非同步观测(卫星轨道、飞机航线)，避免时间插值误差
\item \textbf{误差传播}：通过伴随模式反向传播观测信息至初始时刻
\end{itemize}

{\color{red}\textbf{[伴随模式的代价]：伴随模式需对预报模式每行代码求导，包括物理参数化方案、边界层、辐射等复杂过程。任何模式更新都需同步修改伴随代码，维护成本极高。欧洲中期天气预报中心(ECMWF)花费数十人年开发 IFS 模式的伴随系统。这促使研究者寻找``免伴随''方法，如下节介绍的 DRP-4DVar。}}

\section{降维投影 4D-Var (DRP-4DVar)}

\subsection{DRP 框架避免伴随模式构建}
降维投影(Dimension-Reduced Projection, DRP)方法通过将高维状态空间投影到低维集合子空间，避免构建伴随模式。

准备一组模式空间集合扰动样本 $\mathbf{E}_\mathbf{x} = [\delta\mathbf{x}_1, \delta\mathbf{x}_2, \cdots, \delta\mathbf{x}_K]$，通过非线性模式积分计算得到观测空间集合扰动样本 $\mathbf{E}_\mathbf{y} = [\delta\mathbf{y}_1, \delta\mathbf{y}_2, \cdots, \delta\mathbf{y}_K]$，样本之间近似线性独立。

{\color{red}\textbf{[DRP 的核心思想]：将无限维状态空间投影到 $K$ 维集合张成的子空间($K \sim 50$-$100$，远小于状态维度 $n \sim 10^8$)。假设真实分析增量可由集合线性组合表示：$\delta\mathbf{x} \approx \sum_{k=1}^K s_k \delta\mathbf{x}_k$。这将优化问题从 $n$ 维降至 $K$ 维，避免高维伴随计算。代价是依赖集合质量：若集合未能捕捉主导误差模态，投影会损失信息。}}

将分析增量投影到基向量构成的子空间上，定义 $K$ 维控制变量 $\mathbf{S}$：
\begin{equation}
\left\{
\begin{array}{c}
\delta\mathbf{x} = \mathbf{E}_\mathbf{x} \mathbf{S} \\
\mathrm{H}_i\delta\mathbf{x}(t_i) = \mathbf{E}_\mathbf{y}(t_i) \mathbf{S}
\end{array}
\right.
\end{equation}

通过上式变换得到新代价函数：
\begin{equation}
J(\mathbf{S}) = \frac{1}{2}\mathbf{S}^T \mathbf{E}_\mathbf{x}^T \mathbf{B}^{-1} \mathbf{E}_\mathbf{x} \mathbf{S} + \frac{1}{2}\sum_{i=0}^{n}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{S} - \mathbf{d}_i\right]^T \mathbf{R}_i^{-1}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{S} - \mathbf{d}_i\right]
\end{equation}

在地球系统中，物理状态的维度非常高($n \sim 10^8$)。直接显式构造如此巨大的矩阵在计算上是不现实的。传统方法采用分解技术对 $\mathbf{B}$ 进行简化：
\begin{equation}
\mathbf{B} = \mathbf{U}\mathbf{U}^T,
\end{equation}
然后进行变量转化：$\mathbf{S} = \mathbf{E}_\mathbf{x}^{-1}\mathbf{U}\boldsymbol{\alpha}$，得到
\begin{equation}
J(\boldsymbol{\alpha}) = \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{\alpha} + \frac{1}{2}\sum_{i=0}^{n}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{E}_\mathbf{x}^{-1}\mathbf{U}\boldsymbol{\alpha} - \mathbf{d}_i\right]^T \mathbf{R}_i^{-1}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{E}_\mathbf{x}^{-1}\mathbf{U}\boldsymbol{\alpha} - \mathbf{d}_i\right]
\end{equation}

{\color{red}\textbf{[预条件化的双重作用]：控制变量变换 $\mathbf{S} = \mathbf{E}_\mathbf{x}^{-1}\mathbf{U}\boldsymbol{\alpha}$ 实现两个目标：(1)将先验项简化为 $\frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{\alpha}$，Hessian 矩阵在先验部分为单位矩阵；(2)使 $\boldsymbol{\alpha}$ 在先验分布下标准正交，改善优化问题条件数。这与 3D-Var 的 $\chi$ 空间变换本质相同，但这里通过集合子空间实现。}}

\subsection{兼容非高斯背景场分布假设的 DRP 框架}
传统算法中，$\mathbf{U}$ 是线性变换因此只能表征高斯分布特征。为了表征非高斯特征，可以将线性的 $\mathbf{U}$ 替换成非线性的变换，记作 $\mathscr{D}$。

借助测度论和概率论的相关知识，可以推导出当 $\mathscr{D}$ 为非线性时，同化变分函数会多出一个雅可比行列式项：
\begin{equation}
J(\boldsymbol{\alpha}) = \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{\alpha} + \frac{1}{2}\det\left(\frac{\partial\mathscr{D}(\boldsymbol{\alpha})}{\partial\boldsymbol{\alpha}}\right)^T\left(\frac{\partial\mathscr{D}(\boldsymbol{\alpha})}{\partial\boldsymbol{\alpha}}\right) + \frac{1}{2}\sum_{i=0}^{n}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{E}_\mathbf{x}^{-1}\mathscr{D}(\boldsymbol{\alpha}) - \mathbf{d}_i\right]^T \mathbf{R}_i^{-1}\left[\mathbf{E}_\mathbf{y}(t_i)\mathbf{E}_\mathbf{x}^{-1}\mathscr{D}(\boldsymbol{\alpha}) - \mathbf{d}_i\right]
\end{equation}

{\color{red}\textbf{[非线性变换的雅可比项]：在测度论框架下，从标准高斯空间 $\boldsymbol{\alpha} \sim \mathcal{N}(0, \mathbf{I})$ 到非高斯状态空间 $\mathbf{x}$ 的变换需保持概率测度。变换雅可比行列式 $\det(\partial\mathscr{D}/\partial\boldsymbol{\alpha})$ 修正概率密度，确保 $p(\mathbf{x}) \propto p(\boldsymbol{\alpha}) |\det(\partial\mathscr{D}/\partial\boldsymbol{\alpha})|^{-1}$。代价函数中的 $\det$ 项来源于对数后验的雅可比修正。这为处理偏态、重尾或多峰先验分布提供了理论基础，突破高斯假设限制。}}

\subsection{基于高斯背景场假设的传统 DRP 算法}
利用集合扰动样本估计背景误差协方差：
\begin{equation}
\begin{gathered}
\mathbf{B} = \frac{1}{\sqrt{K-1}}\sum_{k=1}^{K}(\delta\mathbf{x}_k - \overline{\delta\mathbf{x}})(\delta\mathbf{x}_k - \overline{\delta\mathbf{x}})^T \\
= (\mathbf{E}_\mathbf{x}\mathbf{b}_\mathbf{S})(\mathbf{E}_\mathbf{x}\mathbf{b}_\mathbf{S})^T \\
= \mathbf{E}_\mathbf{x}\mathbf{B}_\mathbf{S}\mathbf{E}_\mathbf{x}^T
\end{gathered}
\end{equation}

这里，$\mathbf{B}_\mathbf{S}$ 是 $\mathbf{B}$ 矩阵在由集合扰动样本构成的子空间上的投影，
\begin{equation}
\mathbf{b}_\mathbf{S} = \frac{1}{\sqrt{K-1}}\begin{bmatrix}
1-\frac{1}{K} & -\frac{1}{K} & \cdots & -\frac{1}{K} \\
-\frac{1}{K} & 1-\frac{1}{K} & \cdots & -\frac{1}{K} \\
\vdots & \vdots & \ddots & \vdots \\
-\frac{1}{K} & -\frac{1}{K} & \cdots & 1-\frac{1}{K}
\end{bmatrix}
\end{equation}

{\color{red}\textbf{[$\mathbf{b}_\mathbf{S}$ 矩阵的作用]：$\mathbf{b}_\mathbf{S}$ 是中心化矩阵，将集合扰动 $\mathbf{E}_\mathbf{x}$ 转换为零均值扰动 $\mathbf{P}_\mathbf{x}$。对角元 $1-1/K$ 减去集合平均贡献，非对角元 $-1/K$ 实现去均值化。因子 $(K-1)^{-1/2}$ 确保样本协方差无偏估计。这个矩阵在 EnKF 理论中同样出现，体现 DRP 与 EnKF 的数学联系。}}

利用 $\mathbf{b}_\mathbf{S}$ 对控制变量 $\mathbf{S}$ 进行预处理变换：$\boldsymbol{\alpha} = \mathbf{b}_\mathbf{S}\mathbf{S}$，其中，$\boldsymbol{\alpha}$ 是预处理的控制变量。

投影矩阵：
\begin{equation}
\left\{
\begin{array}{l}
\mathbf{P}_\mathbf{x} = \frac{1}{\sqrt{K-1}}[\delta\mathbf{x}_1 - \overline{\delta\mathbf{x}}, \delta\mathbf{x}_2 - \overline{\delta\mathbf{x}}, \cdots, \delta\mathbf{x}_K - \overline{\delta\mathbf{x}}] \\
\mathbf{P}_\mathbf{y} = \frac{1}{\sqrt{K-1}}[\delta\mathbf{y}_1 - \overline{\delta\mathbf{y}}, \delta\mathbf{y}_2 - \overline{\delta\mathbf{y}}, \cdots, \delta\mathbf{y}_K - \overline{\delta\mathbf{y}}]
\end{array}
\right.
\end{equation}

这里，$\overline{\delta\mathbf{x}}$ 和 $\overline{\delta\mathbf{y}}$ 表示初值扰动样本和观测空间扰动样本的集合平均。$\delta\mathbf{x}$ 和 $\mathrm{H}_i\delta\mathbf{x}(t_i)$ 可以通过投影矩阵被投影到子空间上：
\begin{equation}
\left\{
\begin{array}{c}
\delta\mathbf{x} = \mathbf{P}_\mathbf{x}\boldsymbol{\alpha} \\
\mathrm{H}_i\delta\mathbf{x}(t_i) = \mathbf{P}_\mathbf{y}(t_i)\boldsymbol{\alpha}
\end{array}
\right.
\end{equation}

{\color{red}\textbf{[DRP 与 EnKF 的深层联系]：投影矩阵 $\mathbf{P}_\mathbf{x}, \mathbf{P}_\mathbf{y}$ 精确对应 EnKF 中的集合异常矩阵。两者都利用集合张成误差子空间，区别在于：EnKF 通过卡尔曼增益 $\mathbf{K} = \mathbf{P}_\mathbf{x}\mathbf{P}_\mathbf{y}^T(\mathbf{P}_\mathbf{y}\mathbf{P}_\mathbf{y}^T + \mathbf{R})^{-1}$ 直接更新；DRP-4DVar 在子空间最小化代价函数。前者是``滤波思想''(一步更新)，后者是``变分思想''(迭代优化)，但数学本质相通。}}

通过上式变换得到新代价函数(注意到由于变换函数是线性的，因此此时雅可比行列式为常数)：
\begin{equation}
J(\boldsymbol{\alpha}) = \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{\alpha} + \frac{1}{2}\sum_{i=0}^{n}\left[\mathbf{P}_\mathbf{y}(t_i)\boldsymbol{\alpha} - \mathbf{d}_i\right]^T \mathbf{R}_i^{-1}\left[\mathbf{P}_\mathbf{y}(t_i)\boldsymbol{\alpha} - \mathbf{d}_i\right]
\end{equation}

梯度为：
\begin{equation}
\left(\frac{\partial J}{\partial\boldsymbol{\alpha}}\right)^T = \boldsymbol{\alpha} + \sum_{i=0}^{n}[\mathbf{P}_\mathbf{y}(t_i)]^T \mathbf{R}_i^{-1}\left[\mathbf{P}_\mathbf{y}(t_i)\boldsymbol{\alpha} - \mathbf{d}_i\right]
\end{equation}

{\color{red}\textbf{[免伴随的梯度计算]：梯度公式不再需要伴随模式！观测空间投影 $\mathbf{P}_\mathbf{y}(t_i)$ 由前向模式积分集合扰动得到，$\mathbf{P}_\mathbf{y}(t_i)^T$ 的转置仅是矩阵代数运算。计算复杂度：前向积分 $K$ 个集合成员 $\sim \mathcal{O}(K \cdot C_{\text{model}})$，远小于开发和运行伴随模式的成本。代价是需要存储各时刻的 $\mathbf{P}_\mathbf{y}(t_i)$ 矩阵($m \times K$，$m$ 为观测数)。}}

最小化式上述代价函数得到最优分析。DRP-4DVar 将高维优化问题($n \sim 10^8$)降至低维($K \sim 50$-$100$)，通过前向集合积分避免伴随模式构建，同时保留 4D-Var 的时间一致性优势。

